First-Principles Constitutive Equation for Suspension Rheology 
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We provide a detailed derivation of a recently developed first-principles approach to calculat- 
ing averages in systems of interacting, spherical Brownian particles under time-dependent flow. 
Although we restrict ourselves to flows which are both homogeneous and incompressible, the time- 
dependence and geometry (e.g. shear, extension) are arbitrary. The approximations formulated 
within mode-coupling theory are particularly suited to dense colloidal suspensions and capture the 
slow relaxation arising from particle interactions and the resulting glass transition to an amorphous 
solid. The delicate interplay between slow structural relaxation and time-dependent external flow 
in colloidal suspensions may thus be studied within a fully tensorial theory. 

PACS numbers: 82.70.Dd, 64.70.Pf, 83.60.Df, 83.10.Gr 



I. INTRODUCTION 

Imposing flow on a colloidal suspension distorts the 
microstructure away from that of the quiescent state 
and induces a nontrivial macroscopic stress response. 
The microscopic dynamics underlying this macroscopic 
behaviour are governed by a combination of potential, 
Brownian and hydrodynamic forces which interact in a 
complicated fashion with the solvent flow field, as de- 
scribed by the Smoluchowski equation [l|, Q . The subtle 
balance between these various physical mechanisms gives 
rise to a rich phenomenology, but also serves to compli- 
cate the formulation of tractable theoretical approaches 
Q . In the present work we will present details of a first- 
principles theory, first outlined in j^, which provides a 
unified description of the mechanical response of colloidal 
liquids and glasses to external flow, albeit in the absence 
of hydrodynamic interactions. 

Colloidal suspensions are of great importance for prac- 
tical applications and the performance of many commer- 
cial products and industrial processes often depends upon 
the rheological nonlinearities which occur when colloidal 
particles are added to a Newtonian solvent Q . For exam- 
ple, the phenomena of shear-thinning, shear-thickening 
and yielding are relevant for the even spreading of paints, 
shock absorption in automobiles and the flow of tooth- 
paste, respectively. In order to control and tune the 
rheology of a suspension to meet the needs of a spe- 
cific application it is thus necessary to have an under- 
standing of how the microscopic interactions between the 
constituents influence the macroscopic response Q. The 
challenge to statistical physics is to identify tractable ap- 
proximation schemes which capture the relevant physics 
while remaining simple enough for concrete calculations 
to be performed. 

Quiescent monodisperse systems of spherical particles 
display an equilibrium phase diagram with colloidal gas, 
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fluid and crystalline phases, analogous to those found in 
simple fluids 0- By introducing a sufficent degree of 
polydispersity the mechanisms leading to crystallization 
can be suppressed and the system may undergo dynami- 
cal arrest to a metastable glass or gel phase. Experimen- 
tally, such arrested states are characterized by very slow 
structural relaxation As the system remains non- 

ergodic on measureable timescales it appears as a soft 
amorphous solid with a corresponding elastic modulus. 
Many features of the dynamical arrest observed in col- 
loidal experiments are captured by the mode-coupling 
theory (MCT), both for repulsive glasses for which the 
arrest is caused by steric hinderance jllj and attractive 
glasses and gels where interparticle attraction dominates 
the dynamics [ill, (although ageing dynamics are ne- 
glected within standard MCT). 

It is fair to say that the MCT-based theoretical under- 
standing of quiescent glasses is now in a fairly advanced 
state and that the linear response regime is under con- 
trol. However, much less is known about the nonlinear 
response of arrested states and, despite progress, consen- 
sus remains to be achieved regarding the fundamental 
physical mechanisms at work. A key feature of systems 
under flow is the competition between the timescale of 
structural relaxation, which becomes very large close to 
the glass transition, and the timescale defined by the in- 
verse of the fiow rate. This competition is captured in 
an elegant way by generalized MCT treatments [l3 - [l7j 
for which the memory function generating the slow re- 
laxation becomes reduced by shear via the mechanism 
of wavevector advection. Various nonlinear rheological 
phenomena related to shear thinning are found to arise. 
The earliest generalizations of MCT to treat nonlinear 
rheology focused on the special case of shear fiow HH 
(see also [l^ for a detailed account). More recently, the 
theory has been extended, first to treat time-dependent 
shear [3 ^^'^ subsequently to treat arbitrary time- 
dependent flow These final developments have ele- 
vated the theory to the status of a full constitutive model, 
in the sense that nonsteady three-dimensional flows can 
be addressed. Although full numerical solution of the 
equations in three-dimensions has not yet been achieved. 
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a simplified schematic version of the tensorial theory has 
been developed which provides sensible predictions re- 
garding the flow and yielding behaviour of systems with 
densities close to and above the glass transition [l^ . 

We would like to clarify the relationship of the present 
work with our previous publications on this topic. The 
original outline for extending MCT to treat nonlinear 
rheology under steady shear flow was presented in the 
short paper [l^]- Due to the technical difficulty of nu- 
merically solving the closed MCT equations in three di- 
mensions, tractable schematic and semi-schematic mod- 
els aiming to capture the essential physical mechanisms 
were developed in [2^. The approach of [13] was gen- 
eralized to treat time-dependent shear flow and applied 
to the special case of shear step-strain [l^l- During the 
subsequent development of the theory we made signifi- 
cant technical improvements which involved an alterna- 
tive definition of the transient density correlator from 
that employed in [ll [H [lO] • (The improved definition 
is given by (|68p below.) This new definition ensures that 
the initial decay rate (see (|90| below) remains positive 
for all shear rates, which improves the numerical stabil- 
ity of the theory and leads to simpler and more elegant 
expressions. It is this improved formulation which was 
employed in the short paper Q and which will be ex- 
plained in detail in this work. The simple-shear theory 
recently presented in [l6j also employed the improved for- 
malism and is a special case of the results presented here. 
It has already been solved numerically for steady shear in 
two dimensions leading to results for stresses, structure 
and single particle motion in qualitative agreement with 
Brownian dynamics simulations [2ll . [22j . 

In this paper we present full details of the tensorial 
time-dependent theory that was outlined in [3| ■ Our the- 
ory starts from a well defined microscopic starting point 
(the Smoluchowski equation) and leads directly to an ap- 
proximation for the time-dependent nonequilibrium dis- 
tribution function from which average quantities can be 
calculated. This includes, but is not limited to, an ap- 
proximation for the macroscopic stress tensor as a func- 
tional of the velocity gradient tensor. 

The paper is structured as follows: In Section |ll] we de- 
fine the microscopic dynamics of the colloidal particles. 
In Section IIIII considerations of translational invariance 
are employed to identify the deformation measures ap- 
propriate for describing the system under fiow. In Section 
II VI we develop the integration-through-transients formal- 
ism which generates a generalized Green-Kubo relation 
for the macroscopic stress tensor, the mode-coupling ap- 
proximation of which is outlined in Section |Vl Our ap- 
proximate expression for the stress requires knowledge of 
the transient density correlator. In Section I VI I we de- 
velop a mode-coupling equation of motion for this quan- 
tity which captures the competition between slow struc- 
tural relaxation and the fiuidizing effect of fiow. In 
Section IVIII wc consider application of the integration- 
through-transients formalism to calculation of the dis- 
torted structure factor and show how an integration over 



the anisotropy, which encodes the microstructural dis- 
tortion, recovers the constitutive equation as obtained 
earlier in Section |Vl Finally, in section I Villi we give a 
discussion and outlook to future work. 



II. MICROSCOPIC DYNAMICS 

We consider a system of spherical Brownian parti- 
cles of diameter d dispersed in a solvent with a specified 
time dependent velocity profile v(r,t) — K{t) ■ r, where 
the velocity gradient K{t) is taken to be a fixed input 
quantity. We assume incompressible flow, Tr K{t) = 0, 
and will from the outset neglect hydrodynamic interac- 
tions. Our choice to neglect hydrodynamics is largely 
motivated by the desire to obtain a computationally 
tractable theory, but may also be physically jusifiable for 
high volume fraction states close to the glass transition, 
where particle motion is slow. Given these assumptions 
the distribution function of particle positions obeys the 
simplified Smoluchowski equation [l[ 



dt 



= 0(t)vl/(t), 



where the Smoluchowski operator is given by 

n{t) = Y,9^■[^^-'P^-l^{t)■V,]. 



(1) 



(2) 



In order to simplify notation wc have set both the thermal 
energy fc^T and bare diffusion coefficient Dq equal to 
unity. The total force acting on particle i due to potential 
interations is given by F.^ = —diUN, where C/jv is the 
total potential energy. 



III. TRANSLATIONAL INVARIANCE 

Employing a spatially constant velocity gradient has 
the consequence that the two-time correlation functions 
are invariant with respect to spatial translation. The re- 
quirement that translational invariance holds for the two- 
time correlations provides a method to identify the rele- 
vant affine deformation measures required for our subse- 
quent analysis. A prerequisite is an understanding of the 
invariance properties of the time-dependent distribution 
function ^'(t). If the system is taken to be in thermody- 
namic equilibrium at some initial time to which, without 
loss of generality, can be taken as the time origin tg = 0, 
the Smoluchowski equation ([1]) may be formally solved 
to obtain the distribution at later times 



*(t) 



f* dsn{s) 



(3) 



where we have introduced a time-ordered exponential 
function (Appendix [X| . In the following we will first use 
the formal solution ([3]) to prove the translational invari- 
ance of \E'(t) and then use this result to analyze the two- 
time correlation functions. 
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A. Distribution function 

The fact that K{t) does not depend upon spatial coor- 
dinates suggests that the translational invarance of the 
equihbrium state will be preserved by the Smoluchowski 
dynamics; only relative particle motion is physically rel- 
evant. However, proof of this expectation for a general 
time-dependent flow is complicated by the fact that the 
Smoluchowski operator is itself not translationally invari- 
ant. Shifting all particle coordinates by a constant vector 
r[ = Ti + a leads to the shifted Smoluchowski operator 

n(T',t) = n{T,t) - a - K'^it) - P 
= n{T,t) + A{t), 

where P = di and we have used 

P • K{t) ■ a = a - • P. 



(4) 



(5) 



The positions of all particles are represented by F = 
{ri, • • • jTn}. Substitution of (|4]) into ([3]) gives the dis- 
tribution function following the shift 



(6) 



where we have used the obvious translational invariance 
of the equilibrium distribution vl/g(r') = *e(r). The next 
step is to write the time-ordered exponential explicitly in 
terms of its series definition 



1+ / dsi{n{si) + A{si)) 



J^'; ds {n{s)+Ais)) _ 



+ dsi ds2{n(si) + A{si)){n{s2) + A{s2)) 



+ ■■■■ (7) 

In order to simplify this expression we define the com- 
mutators 



c(2)(si,S2) = [A{si),n{s2)] 

:;(3)(si,S2,S3) = [c(2)(si,s2),f^(s3)] 



(8) 



Due to the fact A{t) is independent of spatial coordinates, 
simple expressions arc obtained for the commutators 

c'-^Hsi,S2) = a - k'^{si)k'^{s2) -P 
c'^^\si,S2,S3) = -a- k^(si)k^(s2)k^(s3) ■ P 



(9) 



Newton's third law states that the sum of internal forces 
in the system must be zero. It follows that 

-Un 



PVe = = E = 0' 



(10) 



where Z is the canonical partition function. We thus 
have 



A{si)'fe{r)=0 

c(")(si,-.- ,S„)*e(r) -0. 



(11) 

(12) 



For each term in the commutation relations ([8]) may 
be used to bring all factors of ft to the left of the inte- 
grand. Application of the relations ([TT]) and ([T2|) then 
eliminates all terms involving k. We are thus left with 



*(r',0 = 1 



dsin{T, Si) 



(13) 



+ ^ dsi^ ds2 n{r, si)n{T, S2) + ■ • ■ ) *e 

exp+ (^J dsn{T,s)j 



(14) 



which is the desired result. We have thus demonstrated 
that the distibution function is translationally invariant 
for any homogeneous fiow field K{t) (although it will be 
anisotropic, reflecting the symmetry of the imposed ve- 
locity gradient). 



B. Two-time correlation function 

We now consider the properties under uniform trans- 
lation of the two-time correlation functions; namely the 
correlation of two wavevector dependent fluctuations. 
Without loss of generality we may express the wavevector 
dependent fluctuations as 



(15) 



where the adjoint Smoluchowski operator can be ob- 
tained from (HI) by partial integration 



(16) 



In contrast to the distribution function, which evolves 
from time t' to time t according to the propaga- 
tor exp_(_ Jj* ds r2(s), the time evolution of fluctuations 
('observables') is dictated by the adjoint propagator 
exp_ Jj, ds r2^(s), where the operation of taking the ad- 
joint reverses the time ordering (see Eas. (jA9[) and (jA10[) ). 
In the present work we will consider only fluctuations 
for which Xf (F) is translationally invariant. This as- 
sumption holds, for example, in the case of the density 
Pq(i), obtained by setting = 1 in p3|) and the less 
familiar case of wavevector dependent stress fluctuations 
<7a/3(q! {^i}) [H, US, obtained by setting 



du{rij) 



dn 



(17) 
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in ([15]). According to this definition aap{<i, {fi}) has di- 
mensions of energy (recall that we set kgT = 1) and 
thus requires division by a volume to become a true 
stress. One should thus bear in mind that the famil- 
iar macroscopic stress arises from taking the average of 
daMq = 0, {rO)/F. 

The two-time correlation function is formally defined 

by 



9k- 



(18) 



As in the previous subsection, we now shift all particle 
coordinates by a constant vector, r,; Vi + a to obtain 
the shifted adjoint Smoluchowski operator 

n^iT',t) = ot(r,t) + a-K'^(t)-p 

= i7t(r,i) + At(i). (19) 

The distribution function appearing in (|18p is invariant 
with respect to this shift as proven in ([T4)). To make 
further progress we employ the identity (|Bip to factorize 
the ordered exponential 



r\dsnUr,s)+A'' (s) (*,ds nUv,s) 
e_ = e_ X 



(20) 



X exp_ 



f* ds'n^ir.s') 



The argument of the second exponential is a complicated 
'interaction representation' of the operator A'^{s), but 
can be simplified considerably using the nested commu- 
tator expansion (|C2[) derived in Appendix [Cj Defining 
the commutators 

5(2)(si,s2) = [n^s^),A\s2)] (21) 

&(3)(si,S2,S3) ^ [r!t(,si),[Ot(s2),At(s3)]] 



we find that 

b^^Hsi,S2) = -a.K^(s2)K^(si)-P 
b'^'^\si, 82,33) = a-K^(s3)K^(s2)K^(si)-P 



(22) 



6(")(si,--- ,s„) - (-l)("+i'a. []/^^(s„_,).P. 



4=0 



Using (PT|) and in the expansion (jC2[) yields a sim- 
plification of the interaction representation operator 

= a • K^(s) (1 + ^ "^^i 

+ ^ dsi^ ds2 k^(.S2)k^(si) + ••■)• P 

d J* ds'n^is')' 



p 



(23) 



where the last equality follows from (|A8|) and the defini- 
tion of the ordered exponential (|A2p . Substitution into 
([20| yields the useful result 



//,dsnt(r,s)+At(s) _ //,dsOt(r,s) 

C — - G X 



(24) 



X exp a 



Substitution of into noting that Pg^. = « k^j^ 

and using the adjoint relation (jA10|) leads to our final 
result 

Cus.it^t') = e-^('i-''(*^*'))-Q^,Jt,i'), (25) 
where the time-dependent wavevector is given by 



k(t,t') = ke 



ftidsK{s) 



(26) 



The condition of translational invariance thus requires 
that the phase factor vanishes and leads us the conclu- 
sion that fluctuation at wavevector k at earlier time t' is 
correlated with one at q = k(i,t') at later time t, as a 
result of the affine fiow. Equivalently, one can employ the 
inverse transformation (see Appendix |X| to show that an 
earlier fiuctuation at wavevector k = q(<, t') is correlated 
with a later one at wavevector q, where 



q(t,t') ^q e 



(27) 



In the following we will refer to (P7|) as the forward- 
advected wavevector and ([26]) as the reverse-advected 
wavevector. 



C. Rheological tensors and advection 



The ordered exponential identified in (\27^ may appear 
unfamiliar, but is nothing more than the inverse of the 
deformation gradient tensor 



E{t,t') 



drjt) 
dr{t') '' 



(28) 



a standard quantity in elasticity theory [2^ [2g. The 
deformation gradient transforms a vector ('material line') 
at time t' to a new vector at later time t. Thus, for 
spatially homogeneous deformations v{t) is transformed 
from the past to the present via 



v{t)^E{t,t')-v{t'), 



(29) 



where Eap = dra/drp, and from the present to the past 
using the inverse deformation gradient 



v{t')^E-\t,t')-v{t). 



(30) 



Taking the time derivative of ([25)1 and applying the chain 
rule to the right hand side generates a differential equa- 
tion connecting the deformation to the velocity gradient 



d_ 

dt 



E{t,t') = K{t)E{t,t') 



(31) 
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Given the boundary condition E{t, t) = 1 the general 
solution of (PT|) is thus 



with inverse (see Appendix E| 



E-\t,t' 



- J*, ds K,(s) 



(32) 



(33) 



The tensors ([5^ and ([55)) are precisely those identified 
from our consideration of the translational invariance 
of wavevector dependent density fluctuations, namely 
k(t, t') = k • E-^{t, t') and k(t, /) = k ■ E{t, t'). 



IV. INTEGRATION THROUGH TRANSIENTS 

A. Nonequilibrium distribution function 

If we now assume that the system was in thermody- 
namic equilibrium in the infinite past the formal solution 
([3]) becomes 



*(i) 



f dsn(s) 



(34) 



Despite being formally exact, the fact that the time- 
evolution is entirely contained within the distribution 
function proves inconvenient when attempting to make 
a closure approximation. For this reason we seek an al- 
ternative, but equivalent, solution for which the observ- 
ables evolve in time and the distribution function remains 
fixed. The procedure is analogous to the passage from 
Schrodinger to Heisenberg pictures in quantum mechan- 
ics and will lead to a Dyson-equation like representation 
of the transient dynamics [l^] . 

Without loss of generality we can separate both the 
Smoluchowski operator and the distribution function into 
equilibrium and nonequilibrium contributions 



*(i) = *e + '5*(t). 



(35) 
(36) 



where 5r2(i) = — ' ('^(O ' ^-^^d <5^'(i) remains 

to be determined. The equilibrium distribution satisfies 
Oe'J'e = 0, whereas nonequilibrium steady states are de- 
termined by rivf = 0. When acting upon the equilibrium 
distribution the nonequilibrium Smoluchowski operator 
generates a term proportional to the stress tensor 



^{t)^e = 50(<)^'e = (K(t):6-)*e 



(37) 



where A : B = X]a/3 -^apBpa and & is the zero wavevec- 
tor limit of the potential part of the stress tensor, which 
we now write in the following form 



r?Ft 



(38) 



In order to avoid any confusion we note that a simple 
transformation to relative coordinates suffices to demon- 
strate the equivalence of p8|) to the more familiar Kirk- 
wood form (|17p [24| . Thus, despite first impressions, the 
compact expression psp to be employed in our calcu- 
lations does not depend upon the choice of coordinate 
origin. 

For the incompressible fiow under consideration we 
have K{t) : S = Tr/«(t) = 0, such that the Kronecker 
delta in (|38l) is irrelevant at this stage of the calcula- 
tion. Substitution of (|35|) and (p6| into ^ thus yields 
an equation for S'${t) 



d_ 

dt 



(5*(t) Vl{t)5'^{t) + {K{t):a)'^^ 



(39) 



which is a simple first order differential equation with 
inhomogencity (/«(t) : o-)5'e. Assuming an equilibrium 
distribution in the infinite past, (j39|) is solved by 



5*(t) 



dti e 



(K;(tl):(5-)*e 



(40) 



as can be verified by either variation of parameters or 
the method of Greens functions. The full nonequilibrium 
distribution function is given by substitution of (j40p into 
(|36p . The nonequilibrium average of an arbitrary func- 
tion of the particle coordinates /({r^}), which need not 
be a scalar quantity, is thus given by 



(Z)"^ = if) + dV 



• r* dsn(s) 



:*T)*e, (41) 



where the integral over F is a phase-space integral over 
all particle coordinates and (•) denotes an equilibrium av- 
erage. The final result, central to the present theoretical 
development, is obtained by partial integration 



^ /*dsf2t(s) 



(42) 



The time development of the observable / is thus gener- 
ated by the adjoint operator. 

As the test function / is arbitrary, it may be removed 
to provide a formal operator expression for the distribu- 
tion function 



*(i) 



r f 

*e + / dti *eK(ii):o-e_i 

J — oo 



(43) 



This expression is equivalent to (j34[) and implies that 
the function to be averaged should be multiplied by 
from the left and integrated over the particle coordi- 
nates. Equation (|43p is the central result of the inte- 
gration through transients formalism and we will demon- 
strate that this provides a very convenient starting point 
for analysis of the nonequilibrium dynamics of colloidal 
suspensions. We note that analogous expressions for the 
nonequilibrium distribution function have been consid- 
ered by Evans and Morriss 27 1 and Chong and Kim [28| . 
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based on the thcrmostattcd SLLOD equations of motion. 
A fundamental advantage of P5)) over ([M)) is that it en- 
ables nonequilibrium averages to be expressed in terms of 
equilibrium averages, under the sole assumption that the 
system was in equilibrium in the infinite past. This leads 
to a Dyson-like representation of the dynamics, reminis- 
cent of time-dependent perturbation theory in quantum 
mechanics. Within our formalism a range of nonergodic 
quiescent states can be generated from this initial state 
by imposing various flow histories. This includes the typ- 
ical experimental protocol for which strong pre-shear is 
applied in order to erase memory and obtain a repro- 
ducible initial state from which the system age may be 
measured. 



B. Exact projection 

A potential pitfall of applying arises from the exis- 
tence of conservation laws which lead to zero eigenvalues 
of and which could, in principle, prevent convergence 
of the time integral. In the following we will construct 
an exact reformulation of (|42p in which slow fluctuations 
are explicitly projected out, ensuring convergence of the 
integral for all values of q and thus preparing the ground 
for subsequent approximations. 

For the Brownian dynamics under consideration only 
the particle number is conserved, dtPq=o = ^^{t)pq=o ~ 
0; which suggests a possible divergence in the hydrody- 
namic limit when the arbitrary test function in (j42p is 
chosen to be a density fluctuation, f ~ Pq- Fortunately, 
this poses no difficulty, as the average in the integrand 
vanishes for all wavevectors 



(44) 



which ensures that density fluctuations do not couple lin- 
early to the nonequilibrium part of the distribution func- 
tion. The fact that ([^^ holds at (7 = follows directly 
from inserting Pq=o = ^ 



(K(t):<Te_i N) 



N{K{t):a) = 0, 



(45) 



where the second equality relies on the condition 
TrK(t) = 0. For nonzero q the vanishing of the aver- 
age in pi)) is a consequence of translational invariance. 
Shifting all coordinates by a constant vector = -I- a 
and applying ([^ leads to 



Jl rfsot(r'.s) 
{K.{t):aeJ Pq) 

-!l,dsn{s) 



(46) 



exp ^— iq • e_ j {K,{t) : cr e''_' pq). 

To prevent violation of translational invariance the aver 
age must be identically zero. Eq. (|44p is thus proven. 

We now introduce the density projector P and com- 
plement Q 

1 



k 



Q = i-P, 



(47) 



(/ql.9q)= dT^J^iT)gqiT). 



where 5*^ is the static structure factor Sq = (p* p^) /N. 
In (|47p we have introduced a bra-ket notation for wavec- 
tor dependent fluctuations. The scalar product of two 
arbitrary variables is defined by an equilibrium average 

(48) 

The projectors (|47|) are both idempotent and orthogonal 
P^=P, Q^ = Q, gP = 0, (49) 
and by construction satisfy 

P\Pci)^\Pci) Q\Pq)^0. (50) 



A direct consequence of (|44| is that our central result 
may be exactly rewritten as 



r' dsnUs) 
}e'l- "Qf), (51) 



(/)'" = (/)+ / dh{K{h): 



where we have again used {K{t) : &) =0. To proceed 
further we consider a time-dependent generalization of a 
standard operator identity (see e.g. [231 the standard 
result) 



j^^dsnHs) _ i\dsQnHs) 



(52) 



+ 



which can easily be verified by differentiation (applying 
the rules from Appendix [A]) . Substitution of this identity 
into ([5T|) and using (jUJ eliminates the complicated sec- 
ond term in (|52p . Finally, exploiting the idempotency of 
Q enables us to write 



f^\dsQnHs) 



e- Q, 



(53) 



which leads to our final result 



J —00 

The projected expression ([5^ is formally equivalent to 
([^^ . The advantage of this reformulation is that the ab- 
sence of linear coupling to density fluctuations is explicit 
in the projected time evolution operator, which ensures 
that approximations to the equilibrium average in (|42p 
will not introduce spurious linear couplings. It should be 
noted that these considerations rely on the incompress- 
ibility condition Tr /« = 0, the relaxation of which would 
indeed lead to a flnite linear coupling for some compress- 
ible flows. 



V. CONSTITUTIVE EQUATION 

A. Generalized Green-Kubo formula 

In order to address the rheology of colloidal suspen- 
sions we now consider application of (|43p to calculation 
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of the macroscopic stress tensor. Replacing / in (|42)) by 
the fluctuating stress tensor elements ([55]) yields a for- 
mally exact constitutive equation 



1 /■* 

a{t) = — dti {K{ti):cTe 

* J — OO 



(55) 



The stress tensor is thus a nonlinear functional of the 
velocity gradient tensor as K{t) appears in the adjoint 
Smoluchowski operator. Equation ((55)) thus goes beyond 
the standard Green-Kubo formulae of linear response 
theory and makes possible a theoretical study of non- 
linear rheology. 

The general result (|55p shows that for arbitrary time- 
dependent flow the integrand must depend upon two time 
arguments, rather than a simple time difference. To illus- 
trate this point, consider the special case of simple shear 
flow (Kapit) = '^{t)5xa5yp, where 7(i) is the shear rate). 
The shear stress is given by 



1 // dsU^is) 

a-xy{t) = / dtij{ti)[—{&xye_^ c 



— OO 

t 



(56) 



where the second equality serves to define the nonlin- 
ear shear modulus, which is a functional of the veloc- 
ity gradient tensor. The lack of time translational in- 
variance demonstrated by G{t,t';[K]) is a general fea- 
ture of two-time correlation functions in systems under 
time-dependent shear and is unavoidable if one wishes to 
go beyond linear response js^l- Time translational in- 
variance can only be recovered by replacing in ([56]) the 
full evolution operator rt'^{t) with the time-independent 
equilibrium operator. For notational convenience we will 
henceforth employ G{t,t') = G{t,t';[K\), omitting ex- 
plicit reference to the functional dependence of the mod- 
ulus on the velocity gradient. 



B. Mode-coupling approximation 

The exact generalized Green-Kubo expression for the 
stress tensor ([55]) requires approximation before explicit 
calculations can be performed. We begin by exactly re- 
expressing (|55p in terms of the projected dynamics (|54p 



1 /■* f' dsQn''(s)Q 

which will be approximated by considering the overlap of 
stress fluctuations with the simplest relevant slow fluctu- 
ations. Due to the projector Q in (j57p the lowest nonzero 
order in fluctuation amplitude is taken to be Pj,Pp, where 
k and p are two distinct wavevectors. We thus define a 
projection operator onto density pairs 



P2 = 



El 

k>p 



1 



PkPp 



(58) 



which we assume will capture the dominant slow fiuc- 
tuations. It should be pointed out that this projection 
operator is only approximate, as we have employed the 
Gaussian ansatz 



(PkPpPkPp) ~ {ptPk){P*pPp) = N'^SkSp 



(59) 



in the denominator [3l| • Although the approximate pro- 
jector (|58p is not perfectly idempotent, it has the advan- 
tage that the normalization is expressed in terms of easily 
calculable equilibrium structure factors. Projecting onto 
density pairs yields an approximation to the stress 

1 /■* f* dsQnHs)Q 

cr(t) = - / dti(K(ti):6-QP2e_^ PaQo"), 

*^ J — OO 



'E T/A 



k>p 
k'>p' 



VN\l_t^ Sk-Sj,,SkS, 



{Pk'Pp'd- 



where the vertex functions are given by 
(K(ti):(TQpk'Pp') 



k'p' 

.(2) 



picPp) 

(60) 

(61) 
(62) 



The second vertex is a tensor, whereas the first is a scalar 
because of the contraction with the velocity gradient. A 
straightforward calculation enables the vertices to be ex- 
pressed in terms of the equilibrium static structure factor 



K'p' = l^{tl)-^P T^r-^k'. 



k'p' 



k dk' 



p' 



(63) 
(64) 



where kp indicates a dyadic product. The wavcvector 
restrictions imposed by the Kronecker deltas appearing 
in and reduce the fourfold sum in (pO)) to a 
double sum. We are thus led to consider the four point 
correlator 



{Pk'P-w e- PkP-w) 



(65) 



In the spirit of quiescent mode-coupling theory we ap- 
proximate the unknown four-point density correlator by 
a product of pair correlators and replace the Q-projected 
dynamics by the full dynamics 



{Pk'P'-k'e- PkP-k) ~ 

{Pk'(^- Pk){Pk'C'- Pk) 



(66) 



For quiescent states this approximate step is well estab- 
lished and amounts to assuming that the density fluctu- 
ations are Gaussian random variables [32|. We assume 
here that this approximation remains valid under flow. 

The essential observation at this point of the calcula- 
tion is that the wavevectors k and k' appearing in (|60p 
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are not independent. The static structure factor and, 
consequently, the vertices (|63p and are clearly invari- 
ant with respect to spatial translation. In order for the 
stress (|60|) to remain translationally invariant we must 
therefore impose that the mode-coupling factorization of 
the four-point correlator ()66p also remain invariant. It is 
this requirement which couples the wavevcctors k and k'. 
Shifting all particle coordinates by a constant vector and 
application of (|24|) reveals that the the right hand side 
of (pS]) is translationally invariant only if k' = k(t,t'), 
where the reverse-advected wavevector is given by ([^5]). 
This identification leads to the appearance of k(t, t') at 
several places within ([60]) and thus generates a highly 
nonlinear functional dependence of the stress on the ve- 
locity gradient tensor. 

Substitution of dM]), dM]) and ^ into using the 
identification k' ~ k(t, t') and replacing the remaining 
discrete wavevector sum by an integral yields 



cr(i) = / dt 



167r3 



kk{t,t') 



X kk- 



ql CI 



Si 



Kit.t')M^ (67) 



where we have defined the transient density correlator 

'i'k(t,t') = ^(Pke^''-^'''^^Vk(M')). (68) 

which remains to be determined. Note that the forward- 
advected wavevector appears naturally in (|68p as a result 
of 



/,UsOt(.s) 

*k(t,t')l*'^ ) - - 



NSk 



(t,t') 



(69) 



The definition in ([BS]) will prove useful later on (in 
Ea.lfTQ])'). when the advection of the wavevector is in- 
cluded into a single time evolution operator. The in- 
fluence of external flow enters (p7)) both explicity, via 
the velocity gradient tensor in the prcfactor, and implic- 
itly through the reverse-advected wavevector. The tran- 
sient correlator captures the slow structural relaxation of 
dense suspensions, thus rendering ()67|) a useful tool for 
the study of arrested states. 

Before proceeding to develop a theory for the transient 
correlator, it is useful to first rearrange (|67| into an alter- 
native form. To this end we introduce the Finger tensor 
[33I , a standard deformation measure from elasticity the- 
ory 



B{t,t') = E{t,t')E^{t,t'), 



(70) 



where the deformation tensor is given by ([32]). The Fin- 
ger tensor describes the stretching (but not rotation) of 
vectors embedded in the material. If we consider the 
quadratic form k-S(t,t')-k and take a derivative with 
respect to i', we obtain 

(k-B(i,t')-k) = \fi{t,t') ■{n{t) + Hi^{t'))-'k(t,t') 



where we have used the definition of the velocity gradient 
pip . Using (|7ip we can express our constitutive equation 
in the following alternative form 



(T{t) = - dt' 




*k(t,t')(*'^')- (72) 



The expression (|72p provided the starting point for the 
development of a simpler 'schematic' constitutive equa- 
tion in 

From ((7^ it can be seen that the stress tensor at any 
given time is built by summing the contributions from all 
possible Fourier modes. Each of these individual contri- 
butions is obtained by integrating up the strain measure 
B{t,t') over the entire deformation history, weighted by 
a decaying, wavevector-dependent memory function. It 
is thus clear that strains occurring much further in the 
past than the longest relaxation time of ^k{t,t') can- 
not contribute to the integral; such distant strain incre- 
ments have been forgotten. Moreover, physical intuition 
suggests that stress contributions at differing wavevector 
should be coupled, as the physical mechanisms driving 
relaxation on different length-scales are not independent 
(a fact ignored in many empirical constitutive equations). 
Within the mode-coupling approximations to be devel- 
oped below this coupling becomes manifest in the nonlin- 
ear functional dependence of the memory function on the 
transient correlator; the equation of motion for $k(i,i') 
is nonlocal in k-space. 



VI. EQUATION OF MOTION 

In order to close our constitutive theory we require an 
equation of motion for the transient density correlator. 
The fundamental approximations to be applied are very 
much in the spirit of standard quiescent mode-coupling 
theory, however, the presence of external flow necessi- 
tates considerable preparatory manipulation of formal 
expressions, not required in quiescent mode-coupling, be- 
fore the theory arrives in a form convenient for making 
closure approximations. The first step on this path is to 
identify the operators which advect the wavevector of a 
density fluctuation (either forwards or backwards) and 
thus capture the purely affine component of the particle 
motion. 



A. Advection operators 

It is a straightforward exercise to show that a density 
fluctuation at advected wavevector k(t, t') can be gen- 
erated from one at k using the exponential advection 
operator 



2k{t,t') ■ K{t) -kit^t') 



(71) 



-f*,ds 5nt(s) 



Pk) - I Pk{t,t')) > 



(73) 
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where the flow term in the adjoint Smoluchowski opera- 
tor is given by 



(74) 



The advection operator describes the purely afBne time- 
evolution of fluctuations which would occur in the ab- 
sence of both Brownian motion and intcrparticle interac- 
tions. It is also useful to identify an anlogous advection 
operator which acts on a density fluctuation to its left 

{Pl\^^'''^'''' = {pk,,,,\, (75) 
where the operator in the exponent is given by 

Wit) = ^ r, • K'^it) ■ (9, + F,). (76) 

i 

The force term in ([75)) arises as a consequence of 
the equilibrium distribution function in our definition of 
the scalar product (|48l) . The operations inverse to ((73)) 
and ((75)) which generate fluctuations at reverse-advected 
wavevectors are given by 



^,...at(.)|^^^ = |pk(M')). (77) 
as can be determined from the inverse relations ((ASp . 



B. Generalized diffusion kernel 

Using ((75)) the transient density correlator ((55)) can be 
rewritten as 



$k(t,t') 



1 



[plU{t,t')p^), 



(79) 



where we have defined the propagator U{t, t') as a prod- 
uct of two time-evolution operators 



i\ — ft' '^^ (^) - ft' '^^ (^) 



(80) 



In order to better understand the properties of U{t,t') 
it is useful to express it in terms of a single exponential. 
Differentiation of ((M| with respect to t and using ((A3() 
yields 



d_ 
dt 



U{t,t') = U{t,t') 



//, ds 5nt (s) I - f* ds sat (^-j 



(81) 



where the subscript on Hl^^{t,t') indicates that this op- 
erator is non-Hermitian. The formal solution of ((81)) is 
given by 



U{t,t') 



(82) 



Our aim to express the transient correlator ((79[) in a 
form suitable for mode-coupling approximation is some- 
thing which can only be achieved with a certain degree of 
foresight, combined with experience regarding the struc- 
ture of the quiescent mode-coupling theory. This has 
the unavoidable consequence that the motivation behind 
some of our formal manipulations can only be fully ap- 
preciated in retrospect. We will attempt to highlight 
these potentially obscure steps when they occur and to 
motivate them as effectively as possible. 

Within the quiescent mode-coupling theory the Her- 
miticity of the equilibrium Smoluchowski operator leads 
to a positive semi-definite 'initial decay rate' Fk = 
— {pl^nlp-k)/NSk determining the short-time dynamics 
of the transient correlator (see Eq. ((89| below). Unfor- 
tunately, the operator ^ll^^^{t,t') appearing in ((82)) does 
not retain this desirable feature and, if not treated care- 
fully, can lead to a time-dependent initial decay rate 
rk(i, t') which changes sign as a function of accumulated 
strain. The resulting zeros of rk(t,t') cause divergences 
in the memory function m]s^(t, s,t') (see the denomina- 
tor of (|103p ) which destroy the numerical stability of the 
theory. These considerations motivate us to seek a refor- 
mulation of the propagator ((80p for which the Hermitian 
operator 

iiUt,t')^JJ'''^^'^nU-^^'''''''^'\ (83) 

becomes responsible for generating the initial decay rate 
of the theory. The operator ((55)) possesses the ap- 
peahng feature that the average {pl^Q,\{t,t')pi^) has the 
same form as that obtained from the equilibrium adjoint 
Smoluchowski operator, except that wavevectors are re- 
placed by their advected counterparts. Positivity of the 
initial decay rate is thus ensured from the outset. 

The operator fll^{t, t') can be brought into the game by 
defining the time-dependent projection operator P(t,t') 
and complement Q{t,t') 



k 

Q{t,t') = 1-Pit,t'), 



(84) 
(85) 



and employing the following exact operator identity to 
rewrite the propagator 



^^Jj,ds'nl^is',t')^^^^^,^JJds'Gis',t')^ 



(86) 



where the operators G{t,t') and H{t,t') are given by 

g{t,t') = e/J"''''''^'^Q{t,t')nle~^^'"''''''^'^ (87) 



The identity ()86|) may be verified by differentiation and, 
despite its complex structure, is simply a further time- 
dependent generalization of (|5^ employed earlier. Note 
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that the advection operator (5ri1^(.s) which is needed to 
build our desired Hermitian operator {ll^{t,t') enters the 
calculation via the advected density fluctuation on the 
right of the projector (|M|) . 

We next take the derivative of (|86|) with respect 
to t and use the lengthy expression obtained to cal- 
culate matrix elements between density fluctuations 
{pI^\ ■ ■ ■ \P]f) /NSk- This calculation yields an equation 
of motion for the transient density correlator featuring a 
generalized diffusion kernel 



(89) 



+ j^dsM^{t,s,t')'^^{s,t') = Ak(i,t')- 

Because matrix elements are taken with respect to a den- 
sity fluctuation at a single wavevector, translational in- 
variance eliminates the wavector sum appearing in (|84|) . 

The initial decay rate rk(t,i') depends upon the Her- 
mitian operator Vl\(t,t') and can be explicitly evaluated 
in terms of the advected wavevector and static structure 
factor 



k{t,t'' 



(90) 



The initial decay rate takes the same functional form as in 
equilibrium, albeit with the static wavevector replaced by 
the advected one. For noninteracting particles only the 
first two terms in ()89|) survive and we obtain the exact 
solution for the correlator in this case, incorporating the 
phenomena of Taylor dispersion (enhanced diffusion in 
the direction of flow) [s^l • 

The generalized diffusion kernel entering (|89p is given 
by the formal expression 



Mk[t, s,t ) = — 



(91) 



where we have defined several new operators in an effort 
to streamline the notation 



BUt,t') = pinlit.t') 
Mt,t') = g{t,t')p^. 



(92) 
(93) 
(94) 



A striking feature of ([OT]) is that the diffusion kernel is a 
function of three time arguments. We recall that for qui- 
escent systems the memory function depends upon only 
a time difference (sij . whereas for systems under steady 
shear it depends upon two time arguments [l6j . The for- 
mer property is a consequence of the time translational 
invariance of the quiescent state (ageing of arrested states 
being neglected within MCT), whereas the latter prop- 
erty is formal recognition of the existence of an absolute 
reference time, namely the time at which the flow was 
switched on. In practice, however, the residual transient 



may be vanishingly small for the steady states of interest. 
The three-time argument memory function (|9ip may be 
interpreted as describing the decay of memory between 
times s and t, incorporating the coupling to the stress 
which is still relaxing from the strain accumulated be- 
tween times t' and t. 

The final term in the equation of motion for the tran- 
sient correlator (|89)) is the 'remainder' term on the right 
hand side, given by 



1 



Ak(t,i') = j:^{piu{t,t',t')g{t,t')p^). 



(95) 



We know rather little about the properties of this term 
and the fact that it remains finite represents a primary 
drawback of the the present approach. (Below we will 
need to approximate it by zero.) Its existence stems from 
our use of the time-dependent projection operator ((85|) 
which only eliminates linear coupling to the density at 
t = t' , leading to Ak(i, t) = 0, but not at differing values 
of the time arguments. 



C. Generalized friction kernel 

Eq. ((M)) represents an exact equation of motion for the 
transient density correlator in terms of a generalized dif- 
fusion kernel. However, experience with quiescent mode- 
coupling theory has shown that it is difficult to accu- 
rately approximate a diffusion kernel and that far better 
results can be obtained when startin g fr om an equation 
of motion involving a friction kernel [l^. We thus seek 
to convert the generalized diffusion kernel Mi^{t,s,t') to 
a generalized friction kernel m-k,{t, s,t'). Following stan- 
dard mode-coupling arguments, the transformation from 
a diffusion to a friction kernel is achieved by splitting 
the time-evolution operator entering the diffusion kernel, 
U{t,s,t') in the present case, into 'reducible' and 'irre- 
ducible' components. As with our previous formal ma- 
nipulations, foresight is required in making this step. The 
specific choice of how to split the time-evolution operator 
is made strategically, with a view to obtain an equation 
of the form (|102p which, together with ((89|) can then be 
conveniently solved to yield (|106p . the desired equation 
of motion. 

The first step is to take the time derivative of U {t, s, t'). 
Using dMl), this yields 



at 



(96) 



We next define a non-Hermitian, time-dependent projec- 
tion operator 

. {plSl'a{t,t')p^) 



with complement Q{t,t') = 1 — P{t,t'), and use this to 
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split Q{t,t') into two contributions 
^U(t, s, t') = U{t, s, t')g{t, t') {Qit, t') + P[t, t')) (98) 

= c7(t, s, t')g"\t, t') + c7(t, s, t')g''"'{t, t'). 

The irreducible operator is simply defined by Q"''{t, t') = 
Q{t, t')Q{t, t'), nothing more can be done with it, whereas 
the reducible operator can be expressed in terms of pre- 
viously encountered quantities 



1 



-AB^{t,t')\ (99) 



Treating the term U{t,s,t')Q'""^{t,t') appearing in 
as an inhomogeneity leads to the formal solution 

U{t,s,t') = U"'{t,s,t') (100) 

+ / dt"u{t",s,t')g'°'^{t",t')u"'{t,t",t'), 



(101) 



where the irreducible time evolution is given by 



which is the solution of the corresponding homogeneous 
equation. Using ([Ml and ([M)) to take matrix elements of 
PIIO)) . (B*(s',t')| • • • \A^{t,t')), and using ^ thus leads 
to an important equation relating the diffusion and fric- 
tion kernels 



Mk(t,s,i') = -rk(s,t')"^k(i,s,t')rk(i,t') 



(102) 



dt" Afk(i", s, t')mi,{t, t", t')Ti,{t, t') 



The generalized friction kernel is given by 

T^^rrwT^' (103) 

^'5^fc(s,t')rk(s,t')rk(i,i') 

which involves the irreducible dynamics. At this stage 
in the calculation the importance of enforcing a positive 
definite initial decay rate becomes clearer. Alternative 
formulations for which advection can cause V\^{t,t') to 
change sign for some values of t and t' will inevitably lead 
to undesirable singularities in the memory kernel ()103p . 
While such singularities could, in principle, be integrable 
and thus remain physical, they would at the very least 
prove inconvenient for numerical implementations of the 
theory. 

We have now almost arrived at our goal to reformulate 
the equation of motion. The final step is to combine ((89)) 
with (|102p . with the former expressed as an inhonioge- 
neous integral equation 



$k(t,t') 



-1 /■* 



This equation can be viewed as a Volterra equation of 
the second kind jssj with an inhomogeneity given by 

..k(M') = j^(|*k(M')-Ak(M')). (105) 

The solution of this equation follows directly from the 
theory of Volterra equations (see Appendix [D| and yields 
the final form for the equation of motion 



(106) 



+ dsmk(t,s,t')^*k(s,t')) =Ak(t,i'), 

where the modified remainder term Ak(t,i') is given in 
terms of previously defined quantities 



Ak(t,t') = Ak(t,t') + rk(i,t') / dsTOk(t,s,i')Ak(s,t')- 



(107) 



Although the integral expression (|107p appears to be 
rather intractable, it can be simplified to a single equi- 
librium average. Substituting the formal solution (jlOOp 
into our expression for the remainder term (j95p recovers 
(|107p with the modified remainder identified as 



Ak(t,t') 



NSk 



{piu'^m,t',t')g{t,t')p^), (108) 



by inspection. The time evolution of Ak(t,t') appearing 
in the equation of motion for the correlator (|106p is thus 
shown to be determined by the irreducible dynamics. 



D. Closure approximation 



(104) 



The equation of motion (|106p we have obtained is for- 
mally exact, but still contains two unknown quantities, 
the generalized friction kernel mk(t, s, t') and the remain- 
der term Ak(i,i'), both of which require approximation 
to arrive at a closed theory. This final stage in our devel- 
opment of a constitutive theory consists of two distinct 
steps. The first of these is specific to systems under ex- 
ternal flow and involves the neglect of certain 'strain en- 
ergy' terms, which, it is hoped, are of minor importance 
in determining the relaxation of the transient correlator. 
Our choice to neglect these terms, which do not arise 
in quiescent mode-coupling theory, represents the only 
point in our development at which an approximation is 
made which lies beyond the established canon of mode- 
coupling projections and Gaussian factorizations. Once 
accepted, this new approximation yields directly an ap- 
pealingly symmetric form for the friction kernel (jl03p 
and causes the remainder term (jl08p to vanish identi- 
cally. The second and final step is then a rather stan- 
dard projection of the friction kernel onto density pairs 
and factorization of the resulting four-point average. 
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We consider first the remainder term (jlOSp . As de- 
scribed in detail in Appendix |El making the assumption 
that the two strain-energy terms introduced there 



E(t,t')= dse 



K {s):cr e_ 



both vanish, gives rise to the simpUfications 



(109) 
(110) 



g{t,t') « gQ{t,t') (111) 

f/--(i,s,i') « U'^^{t,s,t'), (112) 

where QQ{t,t') and U'^^^it, s,t') arc given by (p8)) and 
(|E10[) . respectively We thus find that the remainder 
term vanishes identically 



Ak(i,<') 



NSk 



{plui^^it,t',t')gQ{t,t')p^} = o, (113) 



as both operators in the average are orthogonal to linear 
density fluctuations. 

Our assumption that the strain-energy terms (|109lll0p 
vanish is certainly an uncontrolled approximation, but 
may be partially motivated by considering their depen- 
dence upon strain. For the case of steady flow [l^ , both 
E terms increase in proportion to the strain accumulated 
since startup and thus have an influence which depends 
upon whether the correlator has relaxed completely for 
small strain values. For glassy states our neglect of the E 
terms is thus equivalent to requiring that the yield strain 
remain small, which is consistent with simulation results 
l36l . 113] • In any case, if the yield strain were found 
to be large, then the implicit assumption of a harmonic 
free energy functional underpinnning the entire mode- 
coupling approach would be invalidated, such that the 
new strain-energy terms would be among the least of our 
worries. 

We turn now to the friction kernel p03p . Using the 
previously derived relations dST]), ^ 

and ([M)) we write the numerator in the following, some- 
what more explicit, form 



{Bl{s,t')U'"it,s,t')A^it,t')) = 



(114) 



Approximation (jllOp enables us to replace (517^ (s) in the 
leftmost advection operator by (SOt(s) (see (|E2p ) 



(115) 



which, together with (|112p . can be used to cast the nu- 
merator of the friction kernel, Ea. (|103p . in the more sym- 
metrical form 

^PHs.t')^\ Q(5, t') e- ''''^''^U'^^-{t, t') X 



xe 



//, ds (J^t (s) 



The projected irreducible propagator UQ^{t, s, t') is given 
by ((EToIi and acts in the space perpendicular to linear 
density fluctuations. It is this fact which enables us to in- 
sert, without incurring further approximation, the extra 
projector Q{t^ t') on the left of the irreducible propagator. 
The quantities Q{t, t') pj^.^^ appearing on either side 
of the propagator in (|116p are fluctuating forces which 
do not couple to linear density fluctuations. 

We next introduce the time-dependent projection op- 
erator onto density pairs 



P2(i t') = V l^q(*^*')^p(*.*')^^^^(*.*')^^(*'*')l III 



q>p 



7) 



which is simply (j58p evaluated at the advected wavevec- 
tor and therefore subject to the same Gaussian approxi- 
mation in the denominator. The numerator (jll6p is thus 
approximated by 



X e 



S',dsSnHs)^ ,^ .Mot 



^P2(<,t')Q(i,0f^IPk(M'))' 



(2) 



kqp 



kq'p' 



iV2 



'\Pq{s,t')Pp{s,t') ^ 



X e 



q>p 

q'>p' 

-j;,ds''snf(s') in-, L^dssnHs) 



(118) 



where the vertex functions are given by 



<;,(.M')= 



q(s,t')f^p(s,t' 



NSq[s,t')Sp{s^t') 



(119) 



Employing the factorization approximation for the triplet 
static structure factor -S^qp ~ S'k^'qS'p enables the ver- 
tices to be expressed in terms of known two-point static 
correlations 



^kqp(-^' t') = k(s, t') ■ (q(s, t')c,-(.,t') + P(s, i')cp(.,t')) X 



k,q+p 



(121) 



Kq'p' (^, t') = k'(i, t') ■ (q'(t, t')c,'(M') + P'(i> ^')cp'(M')) X 



XP'5k',q'+p' 



(122) 
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where p = N/V is the particle number density and 
Cq = (1 — l/Sq)/p is the Ornstein-Zernike direct corre- 
lation function. The four-point correlator entering (|118|) 
is approximated in the spirit of quiescent mode-coupling 
theory by factorizing and replacing the Q-projected irre- 
ducible dynamics by the full dynamics 

(123) 

where the Kronecker delta functions enforce translational 
invariance. Using (|118p together with (|12ip . (|122p and 
(|123p to approximate the numerator in (|103p and employ- 
ing the explicit expression (PU)) for the initial decay rate 
in the denominator thus yields our final mode-coupling 
approximation to the memory function 

t') = ^ A ga(M')W)^P(V') ^ (124) 

' ' ^ leW ^ k^{s,t')k^{t,t') ^ ' 

X Vkqp(s, t') Vkqp (t, t')*q(s.t') > s)*p(.s,t') ■ 

The wavevector restrictions enforced by the Kronecker 
delta functions appearing in both the vertices and the 
mode-coupling factorization of the four-point correlation 
function reduce the four-fold sum in (jllSp to a single sum 
over q (which in (|124p has been replaced by an integral) 
and leads to a coupling of the two 'internal' wavevectors 
q and p via the condition p = q — k. 

We have thus arrived at a closed constitutive equation 
for the rheology of colloidal suspensions which requires 
only the number density p, the static structure factor Sq 
and velocity gradient tensor nit) as input. In order to 
calculate the stress tensor one must first solve the non- 
linear integro-differential equation (jl06p for the transient 
density correlator using the approximation Ak(t,i') = 
and the approximate memory function (jl24p with vertex 
functions (|12ip and (|122p . The choice of fiow enters via 
the advected wavevectors. The transient correlator thus 
obtained is then substituted into the constitutive equa- 
tion ([7^ and integrated. There are no adjustable pa- 
rameters in our theory as such, although we accept that 
the approximations we have made are not the only ones 
possible. (In that sense, our theory is one of a larger 
family that could be specified by adjusting parameters 
that remain, at this stage, unidentified.) In the absence 
of fiow the theory for the transient correlator reduces to 
the familiar quiescent mode-coupling theory f3ll |. 

E. Linear response 

It is informative to consider the slow-fiow limit of the 
constitutive equation (j72p for which the stress response 
becomes linear in the velocity gradient tensor. To obtain 



the linear response expression we first neglect the strain 
dependence of the transient correlator and replace ad- 
vected wavevectors by their non-advected counterparts, 
wherever these appear explicitly. We thus obtain 

(125) 

where the quiescent correlator depends only on a time 
diff'erence. We next expand the Finger tensor to linear 
order in the velocity gradient tensor 

B{t,t') = Jl"i^<^)J}ds.-is) 

K 1 + 2 dsK{s), (126) 

where k = {k + k^)/2 is the symmetrized velocity gra- 
dient tensor. Substitution of (|126p into p25p thus yields 
the desired linear response result 

The required quiescent correlator is obtained from solu- 
tion of the K{t) — ?• limit of the equation of motion (|106p 

^$fe(t) + r;,(^$fe(i) + ^ dsmfc(t-s)^$fc(s)^ =0, 

with initial decay rate = k'^ / Sk and memory function 

mu{t) = pJdq^^V^^^^m.it), (128) 

where the vertex function is given by Vkqp = k • (qcg -|- 
pcp). 

The anisotropic part of the integrand appearing in 
(|127p is contained within the factor {•}. The angular 
k-integrals can therefore easily be performed for a given 
K{t) and it follows that for a given steady fiow rate 7 the 
Trouton ratio of extensional viscosity rje = (o'xx — '^vv)/i 
to shear viscosity rjs = (Jxyli takes the values 4 and 3 
when the extensional flow is planar and uniaxial, respec- 
tively. These simple geometrical ratios are in compliance 
with Trouton's rules, familiar to continuum rheologists. 

The elastic limit can be accessed by partial integration 
of (|127p and yields Hooke's law for an incompressible 
isotropic elastic body 

cT{t) = 2G^e{t). (129) 

The infinitesimal accumulated strain is given by 

e{t) = I dsK{s), (130) 

and the single clastic constant predicted by the theory is 
given by the standard mode-coupling expression for the 
plateau value of the shear modulus [2[ 
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For arrested states we recall that the long time limit of 
the correlator $fc(t — > oo) remains finite and serves as an 
order parameter for the transition from a fluid to a glass 
or gel. We note that the present theory makes no predic- 
tion for the bulk modulus, as the assumed compressible 
flow satisfies Tre = 0. 



F. Material objectivity 

The 'principle of material objectivity' expresses the 
requirement that the constitutive relationship between 
stress and strain tensors should be invariant with respect 
to rotation of either the material body or the observer, 
thus preventing an unphysical dependence of the stress 
upon the state of rotation. Material objectivity is an 
approximate symmetry, based on the neglect of inertial 
effects on the microscopic level (i.e. the infiuence of cen- 
trifugal and coriolis forces on particle trajectories). Nev- 
ertheless, many soft materials display this symmetry to 
an excellent level of approximation. The overdamped 
Smoluchowski dynamics underlying our treatment ex- 
cludes inertial effects from the outset and, providing that 
our approximations preserve this, our set of equations 
([72)1 . ((TUel) and (fT24)) should be material objective. 

Material objectivity can be explicity confirmed by us- 
ing and (P7)) to eliminate the advected wavevectors 
in favor of the deformation tensors E{t,t') and B{t,t'). 
When the system is subject to a time-dependent rota- 
tion R{t) the deformation gradient and Finger tensors 
transform as [25| 

E{t,t') = R{t)E{t,t')R^{t'), (132) 
B{t,t') = R{t)B{t,t')R^{t), (133) 

where the hat denotes a tensor in the rotated frame. Ma- 
terial objectivity is verified if the rotated stress tensor is 
found to be given by 

&{t) = R{t)<T{t)R^{t). (134) 

A rather tedious but straightforward calculation shows 
that insertion of the rotated deformation tensors into our 
constitutive equation indeed leads to the relation (|134p . 

VII. DISTORTED STRUCTURE FACTOR 

In order to calculate the distorted structure factor 
'S'k(i) we use the general integration-through-transients 
formula with Q-projectcd dynamics to calculate the 
average of a normalized product of two density fiuctua- 
tions 

S^{t)^^{plp^) (135) 
+ — dt {K{t):(TQeJ QPkPk)- 

J —oo 



In order to arrive at a closed expression for 5*1; we ap- 
proximate the average in the integrand using the time- 
dependent projection operator (jll7|) 

{n{t'):aQjJ''''''''^'^''Qpip^) (136) 

The vertex function appearing to the left of the propa- 
gator is given by ((6T|). albeit with advected wavevectors 
replacing static ones, whereas the vertex appearing on 
the right consists of two terms. The first of these terms is 
identical to (p^ . again with advected wavevectors replac- 
ing their static counterparts, whereas the second term 
is an isotropic contribution which we choose to neglect. 
Our choice to ignore this term is primarily motivated by 
the fact that it seems to be quantitatively small in com- 
parison with the other terms [lB|- (We note, however, 
that the relative smallness of the isotropic term has been 
judged on the basis of calculations performed within the 
grand canonical ensemble [s^ which are not fully con- 
sistent with the requirement of fixed N imposed by the 
Smoluchowski equation.) If we nevertheless choose to ig- 
nore this additional isotropic contribution and perform 
a factorization of the dynamical four point correlation 
function, we arrive at the result 

5k(0 = Sk (137) 

^J_f [ k{t,t') j 

Simple application of the chain-rule enables us to then 
further simplify this expression to arrive at the final form 

dt'^^<i>l^,^,,^{t,t') (138) 

This result has the appealing interpretation that in order 
to calculate the nonequilibrium flow-distorted structure 
factor one has simply to integrate the affinely advected 
equilibrium static structure factor over the fiow history, 
weighted by the transient correlator encoding the struc- 
tural relaxation of the system. We should also note that 
P38p makes clear the difference in philosophy between 
our approach and that of Miyazaki et al. |36l|. In (36j 
the distorted stucture factor is employed as an input to 
the theory. In contrast, we input the equilibrium static 
structure factor, which serves as proxy for the potential 
interactions, and generate the structural distortion as an 
output. 

A notable aspect of our result (|138p is that we can make 
a direct connection between the distorted structure fac- 
tor and our constitutive equation (j72p . By inspection, 
we find that substitution of (|138[) into the following ex- 
pression 
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recovers ([72]) exactly. For the case of shear flow (|139|) co- 
incides with a result of Frerickson and Larson for copoly- 
mers [3^ , in some sense reflecting the Gaussian statistics 
underlying both approaches. 

From (|139|) it can be seen that our choice to ignore the 
isotropic term appearing in (|138p corresponds to neglect- 
ing an isotropic contribution to the stress tensor. While 
the uncertainty regarding this term thus prevents us from 
predicting with confidence the system pressure, it has no 
consequence with regards the off diagonal stress tensor 
elements (shear stresses) or the normal stress differences 
which are, after all, the quantities of most direct rheo- 
logical significance. Moreover, for incompressible systems 
the yielding behaviour of arrested states, as described by 
the yield stress surface [l9| . is invariant with respect to 
isotropic pressure. The nonequilibrium pressure becomes 
relevant when addressing flow-induced particle migration 
and shear banding [i^. Finally, we note that the known 
ensemble dependence of the linear-response Green-Kubo 
result for the bulk viscosity of compressible fluids HO] 
provides circumstantial evidence supporting our choice 
to suppress the isotropic contribution to (|138|) . It seems 
likely that the choice of ensemble will be important when 
employing Green-Kubo-type formulae to calculate diag- 
onal elements of the stress tensor. 



VIII. DISCUSSION 

We have presented a detailed derivation of our recently 
proposed constitutive equation for the rheology of dense 
colloidal suspensions Our theory captures the slow 
structural relaxation arising from potential interactions 
(the 'cage effect') by encoding these in the decay of the 
transient density correlator $k(^)- External flow leads to 
the affinc advection of density fluctuations which com- 
petes with particle caging and tends to accelerate the 
loss of memory. It is this competition of timescales which 
is ultimately responsible for the observed transition be- 
tween a shear thinning viscoelastic fluid and a yield- 
ing glass as a function of the coupling strength. The 
integration-through-transients formalism presented here 
provides a means to calculate stationary averages, cor- 
relation and response functions [igI . |4]| under arbitrary 
time-dependent flow. 

Our consitutive equation (|72p has been simplified in 
to a schematic level and solution of the resulting 
equation for a number of special cases has proved rather 
successful in capturing qualitative aspects of the stress 
response in a robust fashion. So far schematic calcula- 
tions have been performed for steady flows [l^, oscilla- 
tory shear [s^l , step strains [l^, |42| and superposed shear 
and extensional flows [1^ . Of particular note is that the 
schematic constitutive equation predicts a dynamic yield 
stress surface of von Mises form [i^l together with small 
corrections related to the flrst normal stress difference. 
If the schematic model is indeed a faithful representative 
of the phenomenology presented by the full microscopic 



expression (|72|) . as currently seems to be the case, then 
these flndings suggest that ([7^ provides a possible route 
to calculating the yield stress surface from first-principles 
theory. In particular, the details of the yield surface could 
then be related directly to the potential interactions be- 
tween the particles. 

A clear drawback of the present approach is the neglect 
of hydrodynamic interactions between the colloids. To a 
certain extent this may be defensible within the range 
of flow rates for which the basic assumptions of the the- 
ory remain valid, namely Peclet numbers Pe = '^(P / Dq 
less than around unity. Within this regime a simple 
modification of the bare diffusion coefficient may restore 
much of the quantitative error displayed by our Brownian 
dynamics-based expressions when attempting to compare 
with experiment. However, such superficial incorporation 
of hydrodynamics will not affect the yield stress values 
predicted by the theory (occurring in the limit 7 — >■ 0) 
and may thus provide an incomplete picture. A more 
comprehensive method by which hydrodynamics can be 
included into mode-coupling-ty pe t heories is thus desir- 
able (perhaps along the lines of j45ll46j). but care should 
be taken that modifications of the theory do not destroy 
the successful description of the hard-sphere glass tran- 
sition. 

Possibly the most questionable aspect of the present 
closure approximation is our reliance on the approxima- 
tion S = S = (jlOQIllOp . which constitutes a neglect 
of certain stress induced couplings, to arrive at a closed 
equation of motion for the transient density correlator. 
Other than the fact that they increase linearly with the 
accumulated strain, we know rather little about these 
fiow induced terms and one can not draw on experience 
from the quiescent theory to assess their importance. The 
central nature of such an uncontrolled approximation 
within the present formulation could serve as motivation 
to seek an alternative formalism in which they do not 
occur (perhaps employing an operator other than (j83p .) 
In any case, if the yield strain were found to be large, 
then the implicit assumption of a harmonic free energy 
functional, which underpins the entire mode-coupling ap- 
proach, would need correction by higher order terms; it 
seems reasonable to hope that Ak(i,t') can be neglected 
on the same basis. Some further arguments for neglect- 
ing K\s_(t,t'), made in the context of steady shear fiows 
but not specific to that case, are made in ^M,- However 
further work to understand the origin of this term and if 
possible to better justify its elimination would certainly 
be desirable. 
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Appendix B: Operator identity 



Appendix A: Time ordered exponentials 

The positively time-ordered exponential function of an 
arbitrary time-dependent operator A{t) is defined by the 
series expansion [ij SI] 



=1+ dsiA{si)+ / dsi ds2 A{si)A{s2) 
Jti Jti Jti 

dsi ds2 ds3A{si)A{s2)A{s3) + --- ■ (Al) 
ti Jti Jti 

The negatively ordered exponential is similarly defined 
by 

r*2j„4|-„-, /■*2 rt2 rsi 

e_' ^'^1+ dsiA{si)+ / dsJ ds2 A{s2)Aisi) 

Jti Jti Jti 

ri2 psi pS3 

+ dsi ds2 ds3A{s3)A{s2)A{si) + -- - . (A2) 
Jti Jti Jti 

By multiplying out the series expansions in (jAip and 
[2|) the following identities can be proven 



± f*2 dsA(s) =F ff^ dsA(s) , , , 

e+'''i e_ '1 ^ = 1 (A3) 

± Sll dsA{s) T Iti' dsA{s) 



1, 



(A4) 



where care must be taken that causality is respected 
within the multiple integrals. These results are consis- 
tent with the expressions ([5^ and for the deforma- 
tion gradient and inverse, respectively. Using the series 
expansions (|Aip and (|A2p the following useful results for 
the derivatives of time-ordered exponential functions are 
easily proven 

= (A5) 



d ( dsA(s)\ _ _( Jl-- dsA{s) 



dti 



A{h) (A6) 



US - -AitMS (A8) 



Finally we note that the adjoint operation satisfies 

^ M = e_i ^ ' (A9) 
^ M ^ e7 ^ ' (AlO) 



The following operator identity has proved useful in an- 
alyzing the translational invariance of the two-time cor- 
relation functions 



J^^ds{A{s)+B{s)) 



exp_ (^j^ ds A{s) 



(Bl) 



exp_ / ds e^ " B{s)e_ 



where A{t) and B{t) are two arbitrary time-dependent 
operators. The proof proceeds by first defining the oper- 
ator 

A^, s (t^ ds(A{s) + B(s)) 

U+{t2,ti) = ' ^ " (B2) 

^ rt2 pSn — l 

= J2 dsi--- I dSn{A{Sn) + B{Sn)) 
n^oJti Jti 

x---{A{si) + B{si)), 



where it is understood that the n = term is unity. The 
derivative is thus given by 

A[/+(t2, ti) = U+it2,h)[A{t2) + B{t2)]. (B3) 
012 

We now define a new operator 

~ , , ds c(tr>.s) 

U{t2,ti) - ei'i ' (B4) 

= ^ / dsi--- ds„c(t2,s„) • ■ •c(t2,Si), 

where c(t2,s) is at present an arbitrary function. The 
derivative of this new operator is given by 

d 

— J7(t2,il) = U{t2,ti)cit2,t2) (B5) 

012 

Pt2 pSn-i 

+ ^ / dsi--- dSn[c'{t2,S„)---c{t2,Si) 
+ c{t2,Sn) c'{t2,Sn-l) ■ ■ ' c(t2, Si) 



+ C(t2, S«) • ■ • c'(t2, S2) C(i2, Si) 
+ c(t2,S„) • • •c'(t2,Sl)] 

Defining 0(^2, s) (up to boundary conditions) in the fol- 
lowing way 

J-c(<2, s) = C{t2, s)A{t2) - A(t2)c(t2, s) (B6) 
OT2 

leads to many convenient cancellations when substituted 
into (|B5p . resulting in 



^U{t2,h) = U{t2,ti)c{t2M) 
012 



(B7) 



+ U{t2,h)A{t2) - A{t2)U{t2M)- 
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Choosing the boundary condition 0(^2,^2) = ^(^2) yields 



— (7(^2, ^i) + Aih)Uit2,h) - Uit2,h)iAih) + B{t2)). 
0L2 



(B8) 



Multiplying (|B8[) on the right with a negatively ordered 
exponential yields 



(B9) 



U{t2M) {A{t2)+B{t2)). 



Comparison of (|B3p with (|B10p allows the identification 

U+{t2,h)^e^_!' "U{t2,h). (BIO) 

Using (jB2p and (|B4[) to write If^ and [/ as exponentials 
thus yields 

f^^ ds{Ais)+B{s)) f,\'dsA{s) St^dsc(t2,s) 

e_ = e_ e_ . (Bll) 

Finally we require an explicit form for c{t2,s). Using 
our chosen boundary condition 0(^2,^2) — -^(^2) enables 
solution of (IB6I) 



c{t2,s)=e, B{s)eJ , (B12) 



as can be verified by substitution of (jB12p into (|B6I) . 
Substitution of (|B12p into (|Blip thus yields the desired 
result (IbTI). 



lemma, which may be useful in simplifying (|B1|) for cer- 
tain special cases. For arbitrary time-dependent opera- 
tors A{t) and B{t), we find the following nested commu- 
tator expansion 

-Jlds'Ais')^^^^Jlds'A(s')^ - fds, [A{S,),B{S)\ 

J S 

+ /dsi^ ds2[A(si),[A(s2),B(s)]] 

dsi I ds2 / ds3[^(si),[A(s2),[A(s3),S(s)]]] 



+ • • • , (C2) 

which is a continuum version of (jCll) . The proof proceeds 
by first defining the operator 

F{t,s) = e-/^''''^^'"^B{s)J]^''^^''\ (C3) 
Differentiation with respect to t leads to 

^E^ = -A{t)F(t,s)+F{t,s)A{t), (C4) 

where we have used (|A5p and (|A7p . Integration of (|C4p 
from s to i yields a recursion relation 



F(i,s) = B(s) - / dsi(A(si)^^(si,s)-F(si,s)A(si)), 

(C5) 

where we have made the identification _F(s, s) = B{t). 
Straightforward iteration of (jCSP directly generates the 
series ([C2|) . 

Appendix D: Volterra equation 



Appendix C: Hadamard lemma 

For non-commuting, time- independent operators X 
and Y there exists a well-known identity (the Hadamard 
lemma) 

+ + ••• , (CI) 

where [X, y] is the commutator. Proof follows from 
defining the function /(s) = e^'^Fe"*'^, Taylor expand- 
ing about s = and then setting s = 1. (|C1|) can be par- 
ticularly useful in cases for which the sequence of nested 
commutators truncates at a low order, or when the infi- 
nite series can be resummed into a closed form. 

When employing the operator identity (|B1[) a structure 
analogous to the left hand side of (jCl[) arises naturally in 
one of the factorized exponential functions. We are thus 
motivated to find a continuous version of the Hadamard 



A Volterra integral equation of the second kind for the 
function $k;(i,^') can be written the following form (ssj 

<^u{t,t') = g^{t,t') + J^ds Ku{t,s,t')<^u{s,t'). (Dl) 

The wavector subscripts are irrelevant in the formal so- 
lution of this temporal integral equation, but are nev- 
ertheless included here to facilitate comparison with the 
equations appearing in the main text. Subject to certain 
reasonable conditions on both the kernel K]^{t,s,t') and 
inhomogeneity g]s^{t,t'), the solution is given by 

$k(i, t') ^ gi,{t, t') - ds Ruit, s, t') 5k(s, t'), (D2) 

where the resolvent kernel R{t,s,t') satisfies the following 
linear integral equation 

i?k(i,s,i') = -K\{t,s,t') (D3) 

-t- / ds' Ki,{s',s,t')Ruit,s',t'). 
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These results can be used to solve directly the inhomo- 
geneous equation (|104p . Comparing (|104p with (|Dip en- 
ables the identification 

K\{t, s, t') = -Mk(t, s, t')lT^{t, t'). (D4) 

Substitution of (ID4[) into (|D3[) and comparison of the re- 
sulting equation with ()102p identifies the resolvent kernel 

i?k(i, s, t') = -mk(t, s, i')rk(s, t')- (D5) 
Substitution of (|D5]) into (|D2|) yields the solution 

Appendix E: Simplifying the friction kernel 

When attempting to simplify the friction kernel enter- 
ing (|106p it is useful to convert advcction operators in- 
volving (551^ (t) to analogous quantities involving 5^1^ {t). 
This can be achieved using the identities 

j^ass^is) ^ ,i;^'^^^(^)(i + s(M')) (El) 

e;^'^'^^^'"'^^' = [l+nt,t'))e-/^"''^''' (E2) 
where we have introduced the two operators 

Y.{t,t)~ dse_^_ " K {s):(TeJ (E3) 

i/\ /"*! -/'ds'5a+(s') X/ \ - f':ds'5iW{s') 

T.{t,t )^J dse^-'' K-^{s):eTe_ ■ (E4) 

In the main text we will argue that both E and E can 
be set to zero in the final approximation. Proof of ()E1|) 



follows from considering the derivative 

and integration from t' to t. The proof of (jE2p is analo- 
gous. 

The identity (jEip is helpful in simplifying the irre- 
ducible time evolution operator given by 

g'"{t,t')^git,t')Qit,t'). (E6) 

Using (jEip to replace the leftmost ordered exponential 
in (j87p yields a natural division into two terms 

g'"it,t') = {gQ{t,t') + g^{t,t'))Q{t,t'), (E7) 

where 

gQ{t,t')^ J} "'^^''^Q{t,t')nl e-/^' '''' (E8) 

g^{t, t') = JI'^'^^'^Y.it, t')Q{t, t') nl e-/^"'''''''^'\m) 

The purpose of this manipulation is that under the as- 
sumption that Ti{t,t') vanishes we obtain g"'^{t,t') w 
gqit, t')Q(t, t'), which is an operator in the space perpen- 
dicular to linear density fluctuations. This consideration 
then leads us to define the propagator 

which generates dynamics in the space perpendicular to 
linear density fluctuations. 
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